Exchange Monte Carlo for Molecular Simulations with Monoelectronic Hamiltonians 
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We introduce a general Monte Carlo scheme for achieving atomistic simulations with monoelec- 
tronic Hamiltonians including the thermalization of both nuclear and electronic degrees of freedom. 
The kinetic Monte Carlo algorithm is used to obtain the exact occupation numbers of the electronic 
levels at canonical equilibrium, and comparison is made with Fermi-Dirac statistics in infinite and 
finite systems. The effects of a nonzero electronic temperature on the thermodynamic properties of 
liquid silver and sodium clusters are presented. 

PACS numbers: 05.10.Ln, 73.22.Dj, 82.60.Qr 
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In the recent years, the physics of materials and com- 
plex systems has undergone awesome development in the 
field of molecular simulation Significant achieve- 

ments include ab initio, Car-Parrinello molecular dynam- 
ics (MD) linear-scaling ||], or progresses in ergodic 
techniques such as exchange Monte Carlo (EMC) 
also known as parallel tempering, or the Wang-Landau 
algorithm Account of electronic structure is often 
made either implicitely with empirical potentials or more 
explicitely through one-electron methods, either in the 
framework of density-functional theory (DFT) or using 
tight-binding (TB) approximations. In the latter cases, 
the electronic ground state is described via integer occu- 
pation numbers. This is usually adequate for nonmetal 
systems. However the electronic states of metals can 
cross the Fermi surface, therefore integer occupations are 
not appropriate for continuous dynamics. Fractional oc- 
cupation numbers must also be introduced for insulators 
or semiconductors, provided that the temperature is high 
enough for the lowest excited states to be populated. 

Building upon a seminal paper by Mermin |^ who ex- 
tended the Hohenberg-Kohn theorem of DFT to nonzero 
electronic temperatures, several authors proposed to 
combine DFT and fractional occupation numbers 
^. Most of these works were devoted to metals, and they 
tried to find suitable forms of the occupation laws for nu- 
merical stability purposes and a better sampling in the 
Brillouin zone. In particular, it was proposed to replace 
the Fermi-Dirac (FD) function with other expressions 

Even when keeping a FD distribution, the best choice 
for electronic temperature was shown not to be neces- 
sarily related with the nuclear vibrational temperature 
I pof . Additionally, and strictly speaking, the FD function 
holds for an infinite system, but should only be consid- 
ered as an approximation when treating a small molecu- 
lar system such as a metal cluster. 

The goal of this Letter is to show how the true canon- 
ical equilibrium of both ionic and electronic degrees of 
freedom can be simulated through MC methods, for small 
or large sizes. A specific interest in choosing MC meth- 
ods over MD simulations is their greater fiexibility and 
wider range of application. They are very convenient for 
discrete systems, and they can be adapted using statis- 



tical biases to accelerate convergence. They also offer 
a straightforward way to sample grand-canonical ensem- 
bles, which is more difficult with MD [0. As seen below, 
the Monte Carlo method is well suited to the problem of 
electronic thermalization, especially for finite systems. 

We consider the general class of materials modelled 
by monoelectronic Hamiltonians. At any given nuclear 
configuration R, the total energy E depends on the oc- 
cupation numbers {n^}. In the Kohn-Sham formalism, 
these numbers appear explicitely in the expression of the 
density, hence in the energy. In TB models, the band 
contribution is the weighted sum of the one-electron en- 
ergies {si}: -^(R) = ^^niei{'K). The first, most simple 
MC algorithm consists in treating the nuclear (R) and 
electronic (N) degrees of freedom on the same footing, 
by performing random moves on the generalized coordi- 
nates Q = (R, N). Here N = {ui} is the set of instan- 
taneous integer occupation states, which evolve during 
the simulation accordingly with the level statistics. In 
the Metropolis scheme, a change from QoM to Qnow due 
to a change in either R or N is accepted with proba- 
bility acc(Qoid Qnow) = min[l, exp(— /3A_B)] where 
A£; = £^(Qnow) - £^(Qoid)- The MC moves involving 
N must keep constant the total number M of occupied 
states. The only moves of this kind that we consider 
are single exchanges between occupied and unoccupied, 
neighboring levels. For a given ionic geometry R, this 
algorithm clearly converges to the electronic canonical 
distribution in the ergodic limit. 

The local moves involving R and N are now imple- 
mented in the framework of generalized ensembles. In our 
case, we use exchange Monte Carlo by performing simul- 
taneous simulations at various temperatures. With some 
probability p, exchange moves between Qi — (Rj,Ni) 
and Qj — (Rj , Nj ) at the inverse temperatures (3i and 
(3j, respectively, are attempted. They are accepted with 
the probability 

acc(Qi ^ Qj) = min[l,exp(A/3A£;)]. (1) 

In this equation, A/? — (3i — [3j and /S.E — E{Qi) — 
E{Qj). With probability 1 —p, the usual local moves are 
attempted for all trajectories. 

We test this algorithm on a model liquid silver, de- 
scribed by a TB Hamiltonian ||ll|. In Fig. |^, we have 
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plotted the average radial density g{r) at T = 5000 K, 
from three different methods, for a system with 108 Ag 
atoms at constant density p — 10.5 10^ kg.m^''. The 
first method only considers nuclear moves, the electrons 
being frozen at — 0. The second one is the previ- 
ous MC method with adjacent trajectories at 2500 K 
and 7500 K, and p = 0.1. Finally, we perform a sin- 
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FIG. 1: Normalized radial density of liquid Ag at density p = 
10.5 10^ kg.m"^ and vibrational temperature T = 5000 K. Te 
is set either to T or to 0. 

gle trajectory MC calculation using only nuclear moves, 
but fractional occupation numbers given by the Fermi- 
Dirac distribution at Te = T = 5000 K. For each 
atomic configuration, the chemical potential which nor- 
malizes Af is found by a Newton-Raphson minimization 
of the error function x^(/i) = [X]fe"'fe^(M) — A/]^ with 
n™{p) = [1 -I- exp(/3(£fc — p,))]^^, starting from the dis- 
tribution N of the previous MC step. In this case, an 
entropic correction {—TS) to the energy is included with 
the usual form given by Mermin |^ : 

S = -fcB^[nfelnnfe + (1 - n/c)ln(l - n^)]. (2) 

k 

All simulations consist of 10^ nuclear MC cycles. In the 
second method, Af electronic moves are attempted for 
each vibrational move. 

The good agreement between the two calculations with 
Te ^ shows that the present MC algorithm is able 
to equilibrate both the electronic and nuclear degrees of 
freedom, for a quite large system whose electronic statis- 
tical distribution can be safely represented by a Fermi- 
Dirac distribution at the same temperature. 

We now turn to finite atomic metal clusters. Finding 
the electronic average occupation numbers for any given 
nuclear geometry is a combinatorial task, which cannot 
be solved exactly in the canonical ensemble, except for 
very few (< 20) levels. In the bulk limit, the grand- 
canonical ensemble is relevant and gives the FD distri- 
bution. For finite, intermediate sizes, the previous MC 
algorithm can solve this problem numerically. However, 



as is well known for discrete spin systems, exchange "flip" 
moves between occupied and unoccupied neighboring 
states can be expected to be mostly rejected at low tem- 
peratures T < AE/k-B, where AE = -Elumo — £^homo is 
the energy gap between the lowest unoccupied and high- 
est occupied orbitals. The convergence can be greatly ac- 
celerated using the kinetic Monte Carlo (KMC) method 
||l2| of Bortz, Kalos and Lebowitz. We use the KMC 
method in conjunction with the local exchange moves, 
starting from the Tg = electronic distribution where 
only the lowest levels are populated. Tests on small clus- 
ters (up to 16 atoms) have shown that about 10^ KMC 
steps are necessary to ensure convergence towards the 
exact statistical population. In Fig. we plot the aver- 
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FIG. 2: Average occupation numbers for the electronic levels 
of Na4o at T = 500 K. The KMC results are compared with 
the FD statistics, for which the best fit is found at 335 K. 

age occupation number at T = 500 K versus level energy 
for the ground state geometry of Na4o described by a 
TB Hamiltonian |l3) . While the general shape is that of 
a Fermi-Dirac type, the actual FD distribution at elec- 
tronic temperature Tg = 500 K does not match the result 
of the KMC calculation, which is best fitted by a FD law 
at effective temperature — 335 K. 

We can combine the KMC algorithm for the electrons 
with the usual MC moves for the nuclei. To save com- 
putational time, the electronic problem is solved period- 
ically, once every M steps. Again, the EMC strategy 
is used to improve global convergence and reduce quasi- 
ergodicity. However, one must be careful when attempt- 
ing exchange moves between trajectories at different tem- 
peratures, because the energy depends explicitely on the 
nuclear coordinates, but also on temperature via the av- 
erage occupation numbers. The same problem would 
hold for any other temperature-dependent potential, such 
as the effective potentials with quantum corrections used 
in liquids theory [ p^ . 

More precisely, the acceptance probability of an ex- 
change between configurations Ri and Rj initially at the 
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respective inverse temperatures Pi and Pj is now 

acc(Rj ^ Rj) min[l, exp{PjAEj + P^AEi)], (3) 

where AEk = E(Rj,Pk) ~ -B(R,, Pk), k = i oi j. When 
the energies are temperature- independent, Eq. (|I]) is re- 
covered. However, in the present case, they include the 
entropic correction of Eq. (||) and should be considered 
as free energies. The canonical equilibrium of both nu- 
clear and electronic coordinates can be simulated using 
this MC method. Since we are dealing with temperature- 
dependent energies, useful analysis techniques such as the 
histograms methods cannot be applied here. In ad- 
dition, corrective terms to the thermodynamic properties 
appear due to this dependence. Unfortunately, because 
the occupation numbers are obtained numerically, their 
derivatives with respect to P are hard to get. We can rea- 
sonably assume that the corresponding effects are small 
at low temperature, since in bulk metals the electronic 
heat capacity is only a small quantity with a weak (lin- 
ear) dependence upon temperature [|l6[ . 

The previously described MC method is used to simu- 
late the solidlike-liquidlike phase change in small sodium 
clusters. Recent theoretical works Q emphasized the 
significant role of geometric (nuclear) effects on the 
caloric curves. However, the possible effects of electrons 
thermalization have not been considered yet. The influ- 
ence of electronic temperature on shape was discussed 
by Yannouleas and Landman ||l^ , but the jellium model 
used by these authors does not provide information about 
phase changes within an atomistic description. From the 
experimental point of view, the results by Habcrland and 
coworkers are still far from a complete understand- 
ing, especially the complex size-dependence of the ther- 
mal properties. 




FIG. 3: Heat capacity per atom of Nag, obtained from ex- 
change Monte Carlo simulations with frozen (solid line) or 
thermalized (symbols) electrons. The average occupation 
numbers (full circles) are compared with the exact values 
(empty squares) . The results of simple Monte Carlo of Ref. ^ 
are also shown for comparison (dashed lines). 

We first compare in Fig. ^ the heat capacity curves 



of Nas obtained without consideration of electron tem- 
perature, nor use of sophisticated sampling method [po[ , 
with the ones obtained with the present Monte Carlo 
algorithms. Here we can compute the exact average oc- 
cupation numbers by solving the combinatorial problem 
for each configuration, or we can employ the numerical 
KMC procedure, both within the EMC scheme. 35 si- 
multaneous trajectories were propagated with 10^ cycles 
each, the occupation numbers being calculated for each 
configuration. For comparison, the curves obtained at 
Te — but with EMC moves are also plotted on Fig. ||. 
The resulting heat capacities show a very good agreement 
between the two present simulations with nonzero elec- 
tronic temperature. Thus the KMC scheme provides a 
good approach to electronic thermalization. The thermo- 
dynamic curves obtained without using exchange MC, or 
assuming frozen electrons, are very similar. Therefore, in 
this case, we can conclude that (i) the electronic temper- 
ature plays only a small role; and (ii) EMC does not bring 
such an improvement over conventional Monte Carlo. 

This situation becomes somewhat different for larger 
sizes. Na2o, Na4o, Na^g and Nacjg are studied using the 
same methods except the one involving the exact calcu- 




200 

T(K) 

FIG. 4: Heat capacities of sodium clusters calculated from 
exchange Monte Carlo simulations with frozen (solid lines) 
or thermalized (circles) electrons. Also shown for comparison 
are the results of simple Monte Carlo of Ref. ^ (dashed lines). 

lation of the fractional numbers. Again, we also com- 
pare the heat capacities with our previous calculations. 
The whole results are shown in Fig. 0. The overall be- 
havior resembles much that of Ref. 120, with a major 
peak in the heat capacity which marks the onset of the 
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solidlike-liquidlike phase change. However, two signifi- 
cant differences can be noted. First, the melting peak ap- 
pears rather clearly, without any strong premelting fea- 
ture (shoulder or peak) at low temperature. This must be 
contrasted with most other theoretical studies 
which emphasized multistep melting in sodium clusters 
described by various models, but is consistent with exper- 
iments Second, the melting temperature indicated 
by the top of the peak is shifted to lower temperatures by 
about 10-30 K depending on size. Again this brings the 
present results closer to experiments, as our previous ones 
with the same tight-binding model were seen to overes- 
timate melting points in these clusters Electronic 
temperature does not have such a large effect, since the 
Monte Carlo results with Te — nearly match the ones 
with Te = T. This is not quite surprising, since the tem- 
perature required for electronic transitions is large, even 
if it decreases as the cluster grows. However, at least two 
reasons can be invoked to explain the differences between 
the present results and the previously published thermo- 
dynamical curves. Exchange Monte Carlo is known to 
be a convenient method in reducing quasi-ergodicity and 
accelerating convergence [Q. Also, by adding the pos- 
sibility of changing electronic surface, barrier crossing is 
further favored. Thus we expect that the MC methods 
developed here are much more efficient than the simple, 
single trajectories simulations. The new results are then 
consistent with a lower melting point ||23|] . 

Extensions of the present algorithms to mean-field mo- 
noelectronic Hamiltonians other than tight-binding is 
possible. A physical limitation is the relevance of the 



calculated excited levels as single-particle states and the 
neglect of many-body electron interactions. Another 
more important limitation is the difficult combination 
with molecular dynamics, due to the non-explicit de- 
pendence of the average occupation numbers on the nu- 
clear structure, at least for finite systems. The use of 
Car-Parrinello dynamics including entropic corrections 
Ip^ requires some practical approximations, such as the 
Fermi-Dirac form for the occupation numbers, or the as- 
sumption that these numbers do not vary much during a 
short time scale. 

The KMC scheme makes the present method naturally 
suited for use with Monte Carlo sampling of the nuclear 
degrees of freedom. Accelerating procedures 1 20 , psf for 



the total energycalculations in tight-binding models can 
also be a valuable improvement for large systems. The 
methods introduced in this Letter have a wide range of 
applications, for both finite and infinite systems, met- 
als, insulators or semiconductor materials. They pro- 
vide a numerically accurate way of calculating the frac- 
tional occupation numbers, and we gave evidences that 
the Fermi-Dirac statistical distribution is not appropriate 
for small clusters. Combined with advanced Monte Carlo 
techniques such as parallel tempering or the multicanon- 
ical ensemble sampling, these methods enable one to in- 
vestigate the equilibrium thermodynamics of large com- 
plex systems which exhibit various kinds of phase tran- 
sitions due to structural isomerization or electronic exci- 
tations. In this respect, clusters close to the insulator- 
metal crossover or having magnetic properties offer good 
candidates for further investigations. 
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